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Abstract 


A general  labelled  tree  structure  is  introduced  for  a class  of  non- 
uniform  two-dimensional  finite  element  meshes.  The  theoretical  basis  of 
the  structure  and  the  fundamental  access  algorithms  on  the  tree  are  pre- 
sented in  a manner  which  lends  itself  to  extensions  to  higher  dimensions. 

For  use  in  finite  element  computations,  the  tree  is  truncated  considerably 
and  then  the  principal,  relevant  algorithms  are  discussed,  including  the 
refinement  of  the  mesh,  the  computation  of  the  elemental  stiffness  matrices, 
and  the  assembly  and  decomposition  of  the  global  stiffness  matrices  based 
on  nested  dissection  techniques.  An  outlook  to  various  possible  extensions 
of  the  structure  is  also  given. 
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1.  Introduction 

CXirrently  an  experimental  software  system  is  under  development  at  the  ' 

Lhiversity  of  Maryland  which  has  the  following  design  properties: 

(i)  The  system  constitutes  an  applications- independent  finite 

element  solver  for  a certain  class  of  two-dimensional,  linear, 
elliptic  boundary  value  problems  defined  by  a weak  mathematical 
formulation. 

(1.1)  (ii)  Adaptive  approaches  are  employed  extensively.  The  a-posteriori 

estimates  developed  in  [1  ]-[  5]  are  used  to  control  the  adaptive  ! 

processes  and  to  provide  a solution  with  near  optimal  error  j : 

within  a prescribed  cost  range.  i i 

j I 

(iii)  In  the  systems  design  advantage  was  taken  of  the  natural  paral-  f ; 

lelism  and  modularity  of  the' finite  element  method. 

The  general  design  of  the  system  has  been  described  in  [ 6 ] . Since  i 

it  represents  an  experimental  prototype  rather  than  a production  system, 

extensive  provisions  for  evaluating  the  performance  are  incorporated.  ! , 

A principal  feature  of  (1.1) (ii)  is  an  adaptive  mesh  refinement  algo-  i i 

• i 

rithm.  Briefly,  after  a solution  has  been  obtained  on  some  mesh  error 

i 

indicators  are  evaluated  on  the  individual  elements  and  from  these  a very  j 

reliable  estimate  of  the  error  in  the  energy  norm  is  composed.  The  error  I j | 

is  (asymptotically) optimal  for  the  degrees  of  freedom  used  if  all  indica-  I • | 

tors  are  essentially  equal.  This  provides  the  basis  for  the  refinement  | | 

algorithm  which  in  essence  divides  certain  elements  so  as  to  achieve  a more 
eqxoal  distribution  of  the  error  indicators  (see,  e.g.,  [5J).  j 


The  efficiency  of  such  a refinement  strategy  depends  critically  on 
the  design  of  the  data  structure  for  the  meshes.  This  is  the  topic  of 
the  present  paper.  More  specifically,  we  present  here  a tree  structure 
for  the  class  of  meshes  considered  in  [6].  Chapter  2 outlines  the  general 
definition  of  the  meshes  and  introduces  the  theoretical  basis  of  the  tree 
structure  for  them  and  the  fundamental  access  algorithms  on  the  tree.  Tliis 
structure  is  not  restricted  to  two-dimensional  meshes,  but  the  extension 
shall  not  be  considered  here.  For  the  finite  element  computations  ' cider 
discussion  the  tree  can  be  truncated  considerably.  This  is  discussed  in 
Chapter  3 together  with  the  principal  algorithms  needed,  namely,  the  refine 
inent  of  the  mesh,  the  confutation  of  the  elanental  stiffness  matrices,  and 
the  assembly  and  deconfosition  of  the  global  stiffness  matrix.  Finally, 
in  Chapter  4 we  indicate  some  possible  extensions  of  the  structure. 
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2.  The  Basic  Data  Structure 
2.1  Domain  aiid  Mesh  Definition 

The  parallelism  of  the  design  property' (1.1) (iii)  is  on  the  procedural 
level  rather  than  the  instructional  level.  It  is  specified  in  terns  of 

i 

* processes  which  are  autonomous  units  with  their  own  programs  and  data. 

These  processes  run  in  parallel  and  communicate  asynchronously  in  a 

I 

limited  and  highly  structured  manner. 

f A natural  parallel  process  structure  for  the  system  derives  from  the 

I familiar  substructure  analysis  in  engineering  design.  The  domain  ffi  c r'" 

t is  defined  as  the  union  52  = ^^^11^0  11  ...  U of  finitely  many  closed, 

I bounded  subsets  52^  c R which  have  nonempty  interiors  52^  such  that 

! 52^  n 2-  = ^,  i j . On  each  subdomain  2.  a finite  element  mesh  is 

! introduced  and,  to  a considerable  extent,  the  confutations  on  the  different 

' subdomains  are  performed  in  parallel  (see  [ 6 ] ) . 

' For  the  design  of  a reasonably  efficient  mesh  refinement  algorithm 

some  restrictions  on  the  choice  of  the  subdomains  are  desirable.  It  is 

assumed  that  each  2^  is  a diffeomorphic  image  of  some  fundamental  figure 
! 2 

, ) F in  R on  which  a sinfle  hierarchy  of  subdivisions  can  be  defined.  On 

» the  intersections  of  the  subdomains  these  diffeomorphisms  have  to  satisfy 

appropriate  compatibility  conditions.  The  finite  element  meshes  on  each 

• ^ 

2^  consist  of  curvilinear  elements  which  are  first  defined  on  F and  then 
mapped  into  2^.  Thus  the  mesh  construction  takes  place  in  F and  for  the 
' discussion  of  the  data  structures  it  suffices  to  restrict  the  attention  to  F. 
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From  a practical  viewpoint  there  are  essentially  only  two  types  of 
fundamental  figures  F that  should  be  considered  here,  namely,  a square 
or  an  equilateral  triangle.  In  order  to  keep  the  data  structures  manage- 
able, it  is  advantageous  to  require  that  the  subdivisions  of  the  chosen 
figure  F consist  solely  of  the  same  type.  For  instance,  if  F is  the 
closure  of  the  open  unit  square 

(2.1)  Qq  = {X  € I 0 < < 1,  i = 1,2}  , 

tlien  admissible  meshes  on  may  be  defined  as  collections  M of  closed 
squares  in  Qq  which  are  generated  by  recursive  application  of  the  tw 
rules; 


(i)  The  mesh  M consisting  only  of  itself  is  admissible 
(ii)  If  M is  an  admissible  mesh  on  Q^,  then  the  mesh  M'  is 
(2.2)  admissible  that  is  obtained  from  M by  subdividing  any  one 
closed  square  Q of  M into  four  congruent  squares  of  half 
the  side  length  of  Q. 

A typical  mesh  on  Qq  generated 
iij  this  way  is  shown  in  Figure  1.  Clearly, 
the  refinanent  introduces  "irregular"  points 
--marked  by  small  circles --which  are  not 
comers  of  all  the  squares  incident  with 
them.  If  conforming  elements  are  used, 
the  solution  at  these  points  is  specified 
by  continuity  conditions;  this  results  in 
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as  carriers.  For  simplicity  of  the  presentation  we  restrict  ourselves  to 
fJermitian  elements  for  which  all  degrees  of  freedom  are  concentrated  in 
the  comers  of  the  squares.  Lagrangian  elements  requiring  additional  nodes 
could  be  used  as  well;  they  increase  only  slightly  the  level  of  complexity 
of  the  algorithms. 

2.2  The  Basic  Tree  Structure 

A widely  used  data  structure  for  finite  element  computations  is 
based  on  a list  of  the  nodes  each  pointing  to  the  elonents  to  which  it 
belongs,  and  of  a list  of  the  elements  which  in  turn  point  to  the  nodes 
incident  with  them.  This  essentially  static  structure  is  not  very'  effi- 
cient when  mesh  refinements  are  introduced.  Instead,  the  recursive  defini- 
tion (2.2)  of  the  admissible  meshes  suggests  the  use  of  a tree  structure 
that  corresponds  to  the  refinement  process  and  has  several  obvious  advan- 
tages. 

Evidently,  the  subdivision  of  a square  in  some  mesh  recjuires  only  a 
simple  extension  of  the  tree  while  in  the  node/element  list  structure 
various  changes  are  needed  in  widely  dispersed  places.  VsTien  the  tree 
becomes  too  large,  it  is  easily  partitioned  into  logically  coherent  parts 
for  storage  on  secondary  devices.  Since  the  tree  structure  reflects  the 
refinanent  process,  it  provides  for  an  efficient  deconqxjsition  of  the  global 
stiffness  matrix  corresponding  to  the  well-knowu  nested  dissection  tech- 
nique (see,  e.g.,  [7  ]).  The  tree  structure  also  allows  for  a rather 
efficient  treatment  of  the  irregular  points  discussed  above  in  connection 
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with  Figure  1.  In  particular,  it  turns  out  that  it  suffices  to  represent 
only  the  regular  points  on  the  tree. 

In  this  section  we  introduce  a fonnal  definition  of  the  basic  tree 

structure  which  will  then  be  simplified  in  Chapter  3. 

12  IT 

Let  B = {e  ,e  } be  the  set  of  the  standard  basic  vectors  e = fl,0) 

T 2 

and  e*"  = (0,1)  of  R . For  any  one  of  the  four  subsets  E c B the  cardinality 
is  denoted  by  |E|.  In  the  |E| -dimensional  plane  u+span  E the  open  "square" 

Q with  center  u and  side  length  h > 0 is  defined  by 

(2.3)  Q = sq(u,h,E)  = {x  ? | ||x-u|l  < 4 h,  x^e^  = u. , V e E}  . 

112 

K’e  write  dim  Q = |E|.  In  particular,  sq(y(e  +e  ),1,B)  is  the  open  unit 

2 

square  (2.1)  and  for  any  u € R and  h > 0,  sq(u,h,0)  is  the  point  u. 

Now  let  A be  the  set  of  nine  vectors 

(2.4)  A = (a  € R^  I a.  = -1,0, +1;  i = 1,2} 
and  set 

(2.5)  Aj  = {a  € A I a.  = 0 if  e^  E},  V E c B , 

(2.6)  E[a]  = {e^  € E I a.  = 0},  V a € A . 

Then  the  faces  of  any  open  square  Q = sq(u,h,E),  E c B,  are  the  3*^™^  squares 

(2.7)  P = sq(u  + 2 ha,  h,  E[a]),  V a € Ag  . 

Note  that  Q itself  is  the  unique  |E| -dimensional  face  of  itself. 
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In  line  with  the  refinement  rule  (2.2)(ii),  we  introduce  a subdivision 
operator  o which  associates  with  any  open  square  Q =*  sq(u,h,r),  I;  c B, 
the  set  a(Q)  consisting  of  the  squares 

(2.8)  sq{u  ha,  T h,  I;  ~ l;[a]),  V a ( A^.  . 

Note  that  o(Q)  contains  exactly  one  jxjint,  namely,  sqfu,  4 h,  0).  for  H = 0 

this  is  the  jioint  Q itself;  othenvise,  it  is  tne  center  of  Q. 

In  order  to  formalize  the  refinement  algorithm  (2.2)  on  the  closed  unit 

1 2 

sqiure  we  mbed  in  the  larger  open  square  Q , = sq(e^+e‘',4,B) . 

Now  we  establish  a fixed  rudimentar)’  tree  as  follows: 

(i)  The  root  of  represents  the  square  Q 

(ii)  The  successors  of  the  root  are  four  nodes  corresjionding  to  the 

(2.9)  squares  sq((e  +e  )+a,2,B'F.[a])  of  o(Q  with  a € A^,  a s 0.  ^ 

(iii)  .Any  node  constructed  in  (ii)  representing,  say,  Q = sq(u,2,E), 

luis  successor  nodes  representing  the  squares  of  a(Q) 

with  a ( Aj;,  a i 0. 

It  is  easily  seen  that  the  resulting  rudimentar>'  tree  Tq  has  exactly 
nine  terminal  nodes  representing  the  open  unit  square  Jind  its  eight  onc- 
;uid  zero- dimensional  faces. 

Now  all  admissible  subdivision  trees  T are  obtained  recursively  by 
application  of  the  two  rules: 

Tl  7 

We  use  here  the  standard  partial  ordering  on  R"  defined  by  x 2 0 when- 
ever Xj^  2 0,  i ” 1,2. 
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(i)  Tq  is  an  admissible  tree. 

(ii)  Let  r be  an  admissible  tree  and  consider  any  terminal 

I node  of  T corresponding  to  a two-dimensional  square  Q. 

(2.10)  Then  a new  admissible  tree  T'  is  obtained  if  we  attach 

to  each  terminal  node  of  T that  represents  a face  P i 

• of  Q exactly  successor  nodes  corresponding  to  | 

j 

the  squares  (2.8)  of  a(P).  i 

I 

For  any  such  tree  T the  terminal  nodes  representing  two-dimensional 
squares  correspond  exactly  to  the  interiors  of  the  undivided,  closed  squares  : 

of  an  admissible  mesh  defined  by  (2.2).  ; 

’ I 

Clearly,  these  admissible  trees  grow  rapidly  very  large;  hence,  for  the 

I 

! imjilementation  they  have  to  be  reduced.  Tnis  will  be  discussed  in  Chapter  3 

I 

' below  where  also  some  illustrative  examples  are  given. 

2.3  Labels 

Let  T be  an  admissible  tree  as  defined  by  (2.10).  We  label  the  nodes 

I 

of  T as  follows: 

(i)  The  root  of  T is  labelled  X(root)  = (1,1)^. 

I , 

(ii)  For  any  nonterminal  node  of  T representing  the  square 

(2.11) 

Q = sq(u,h,E)  each  successor  node  p corresponds  to  a 
square  P = sq(u  + ha,  j h,  B^E[a]),  a € Ag,  of  o(Q). 


This  node  p is  labelled  X(p)  = a. 
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These  labels  allow  for  an  easy  reconstruction  of  the  square  Q = sq(u,h,F.) 
corres{X5nding  to  a node  q of  T,  In  fact  let 

(2.12)  rr  = path  (q.root)  = . . . .p">,  P^  = q.  = root,  n > 1 , 

be  the  unique  path  from  the  node  q to  the  root.  Suppose  that  p 
corresixjnds  to  the  square  Pj^  = sq(u  ^ “ l,...,n.  Then  we  have 

by  construction 

' t n , , n.  1 2 

h =4,  u =\(p)=e  +e 
n 

^k-1  “ i . k = n,n-l , . . . ,2 


and  hence 


n-1 


(2.13)  h,  = u*^  = X(p")  + I X(p^)2^’‘^’^'^\  k = n,n-l,...,l 

^ j=k 

Moreover,  by  (2.8)  we  have 


1 


(2.14)  Ej^  = E ~ E[X(p’^)],  k = n,n-l 1. 

The  labelled  tree  T allows  also  a sin^jle  reconstruction  of  the  faces 

I 

and  neighboring  squares  of  a given  square.  For  this  it  is  useful  to  intro-  j 

duce  some  notation.  I 

1 

Note  first  that  the  set  A of  (2.4)  is  a multiplicative,  commutative 


semigroup  under  the  product 


r 


X-.  ..JHjLiJ^iiPfwilPM.Ji  LI  111  — 
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With  the  members  ot  A we  form  strings  of  finite  length 

(2.16)  a = a[l]a[2] . . .a[n] , a[i]  € A,  i = l,...,n  . 

For  such  strings  the  following  three  operations  are  well  defined; 

substring:  a[i:j]  = a[i]ci[i+l]  • • •a[j  ] , 1 ^ i £ j £ n 

(2.17)  concatenation:  a • p = a[l] . . .a[n]p [1] . . .p fm] 
multiplication:  a o*  a = (a'.->a[l])  (a<ia[2]) . . . (a‘^[n]),  V a t A . 


Now  suppose  that  on  the  subdivision  tree  T the  node  q represents 
the  two-dimensional  square  Q = sq(u,2^  ^,B)  with  side  length  h = 2^ 


n i 3.  We  wish  to  find  the  nodes 
of  T corresponding  to  the  comer 
{K)ints  i,j  = il,  the  sides 

i = ±1,  j = 0,  or  i = 0,  j = ±1, 
and  the  neighboring  squares  of  the 
same  size  S...  i = ±1,  i = 0,  or 
i = 0,  j = ±1,  if  they  exist  (see 
Figure  3.) 


V 

S.-. 

S ’■.c: 

i Q 3 

s 

f.c 

p. 

- • 

9 

J ^ 

Figure 

it  is  obvious  that  the  centers  of  these  various  "squares"  are  given  by 
u(P.j)  = u(Q)  + (J),  i,j  = ±1 


(2.18) 


u(Tjj)  = u(Q)  + 2^'"  (j)  1 i = ±1.  j = 0 

> and 

u(S..)  = u(Q)  + 2^'*^  (j)  J i = 0,  j * il 
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Let  n be  the  path  (2.12)  from  q to  the  root  and  a(q)  = string(tT)  the 

corresponding  label  string  with  a[k]  * X(p*^)  - k = 1 n. 

For  ease  of  notation  we  denote  by  co  any  one  of  the  conponcnts  of  u(Q). 
By  (2.15)  we  then  have 

xl,2-n  . 2_3-n  ,-l.n-2  ,n-l  .n 

u)»X2  ♦X2  ♦...  + 2X  +X  +X  . 


An  addition  or  subtraction  of  2^  evidently  generates  a carry.  >bre 
specifi'cally,  assune  that 


. .k-2  .k-1  .k-2  . ^ , 

X=X»...  *X  ,X  =-X  ,k22, 


then  it  follows  readily  that 


. 1,2-n  . 2 3-n  ^ . .n 

CO-X2  »X2  +...+X 


(2.19) 


♦ X^2^''^  = ^k2k-(n-l)  ^ ^ ^n 


CO  - X^2^’’^  = -X^2^’^  + X^2^'’^  + ...  + X*' 


. .l,3-n_  , lo2-n  ,k-2^k-n-l  ,k-l-k-n 

CO  “Afc  “•••*A  » *A  • 


In  order  to  translate  this  into  label  strings  we  introduce  for  any 
node  q with  label  string  a = string (n)  the  indices 


(2.20)  k.  - k.(q)  - 1 + min(k|X^  - -xj) , i = 1,2,  a[j] 


( 

( 


Since  always  1,  * -1,  ^ ” 1»  i“  1»2,  the  existence  of  these 

indices  is  guaranteed  for  any  node  q for  which  n = path (q, root) 
has  at  least  length  n > 2.  Evidently  the  nodes  with  n = 0,1  arc  nonterminal 
nodes  of  the  rudimentary  tree  and  have  no  direct  bearing  on  our  subdivision 
process. 

Now  let  p..,t..,  and  s..  denote  the  nodes  of  T (if  they  exist) 

‘ ij  ij  ij 

corresponding  to  the  "squares"  , T^.,  and  , respectively.  Moreover, 

T 

assume  that  X(q)  = (Xj,X2)  . Then  we  obtain  from  (2.19)  the  following 
formulas  for  the  label  strings  P(p^j)  = string (path (p^^ , root)) : 

P(P.x  .X  ) - a[2:n] 

P(P.Xj,X2^  = {e^»ar2:k2-l]}  * a[k2:n] 

P(Px  -x  ) “ {e^«)a[2:k,-l] } • a[k.  :n] 

(2.21)  ^1*  ^2  ^ ^ 

|'{eW[kj^:k2-l] } • a[k2:n]  if  kj  < k2 


P(PXi^X2^  = j a[k2:n] 


if  kj  * k^ 


|^{e'^<»a[k2:kj^-l] } • a[kj:n]  if  k^  > k. 
Similarly  it  follows  for  P(t  — ) = string  (path(t^^  ,root))  that 

P*^^0’-X2^  =*  • a(2:n] 

P(t.x  ,0^  “ ^2^'^  ’ “t2:n] 

(2.22) 

Pttp.x^)  * {eW[l;k2-l]}  • a[k2:n] 

P(tx^^0^  " ^e^^^a[l:kJ^-l]}  • o[kj:n] 


i 
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cuid  for  P(s^j)  » string(path(s. j ,root))  that 

■ (a‘)  • 

PfSp.X-,^  * ((_J)>5a[l:k2-l]}  • a[k2:n] 
= l(‘}>-)a[l:kj-l]}  • a[kj:n] 


Figure  4 shows  the  nodes  , t^j , and  on  the  tree  and  the 

paths  represented  by  the  various  label  strings.  In  order  to  find  the 
different  nodes,  we  need  only  follow  the  indicated  paths  in  T given  by 
the  label  sequences  written  besides  than. 

The  following  algorithm  retrieves  the  nodes  and  their  labels  on  the 

path  n = path(q,root)  in  the  arrays  p and  a,  respectively.  It  also 

gives  the  length  n of  n and  the  two  indices  k^  = k^(q),  i = 1,2. 

Algorithm  UfJ-Path 
Input  (q) 
n 1;  p[l]  :=  q 

kj  :=  0;  a[l]j  :=  X(q)^,  for  i = 1,2 
r father (q) 
hTiile  r ^ nil  do 
begin  n :=  ml;  pfnj  :=  r 
For  i =■  1,2  do  begin  a(n]j^  :=  X(r)^; 

if  (k^-0)  and  (a[n]^^X(r)j)  then  k^  :»  n+1  end 
r father (r) 

end 


Output  (n,p,a,k) 
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Here  father (r)  returns  the  father  node  r if  r is  not  the  root,  else  it 
returns  "nil". 

Once  the  arrays  p and  a have  been  established  in  this  way,  the 
following  algorithm  follows  the  pxjrtion  of  the  down-path  starting  at  the 
node  p[k]  = p defined  by  the  label  sequence 

4 {a[k]a[k-l]...ci[<]},  k>  ^ 2 1 . 

It  returns  the  endpoint  of  the  requested  part  of  the  down-path  if  that 
point  exists,  else  it  returns  "nil". 

Algorithm  Down-Path 
Input  rp,a,k,«, 3^,32) 
i :=  k;  r :=  p[i] 

While  (ii^nil)  and  (i>^)  ^ 
begin  i:=  i-1; 

r :=  son(r,3ja[k]j,3  2a[k]2)  end 
Output  (r) 

T 

Here  son(r,Xj^,X2)  returns  the  descendant  of  node  r with  label  (Xj^,X2) 
if  it  exists;  else  it  returns  "nil". 

In  connection  with  the  treatment  of  irregular  nodes,  we  shall  need 
the  smallest  open  squares  in  the  mesh  which  contain  both  the  given 

two-dimensional  square  Q and  its  open  side  T.j,  i “ ±1,  j ■ 0,  or 

i = 0,  j = +1.  If  q is  again  the  node  of  T corresponding  to  Q and 

T 2 

of  our  halving  strategy- - the  father  p 

of  q always  represents  the  squares  Q ^ n * Qn  \ • Nbreover,  it  is  geo- 

metrically  obvious  that  p ^ and  p ^ correspond  to  q and  Qq  ^ , 

respectively.  Thus,  these  four  cases  are  a simple  byproduct  of  our  algorithm. 


! 
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3.  Implementation  Aspects 

3.1  The  Truncated  Subdivision  Tree 

As  indicated  earlier,  for  the  finite  element  confutations  considered 
here  it  is  unnecessary  to  implement  the  full  subdivision  tree  T discussed 
in  Chapter  2.  In  fact,  we  may  truncate  T considerably  by  deleting 
branches  carrying  information  not  needed  in  the  calculations.  This  trunca- 
tion is  achieved  in  several  steps. 

When  an  open  square  Q = sq(u,2h,B)  is  subdivided,  a node  p represent- 
ing the  center  point  u = sq(u,h,|))  is  introduced  in  T.  Thereafter,  if  one 
of  the  squares  incident  with  u is  subdivided,  p receives  one  descendant 
node  representing  again  the  point  u = sq(u,h/2,|)) . This  repeats  itself; 
that  is,  a string  of  nodes  with  single  descendajits  is  generated,  all  of 
which  represent  the  center  u of  Q.  A first  truncation  of  T therefore 
consists  in  deleting  all  nodes  that  represent  a point  u = sq(u,h,|9),  h > 0. 

This  leaves  us  only  with  nodes  on  T that  represent  sqxoares  of 
dimension  one  or  two.  Once  such  a square  has  been  subdivided,  it  no  longer 
plays  a direct  role  in  the  computation.  Hence  from  now  on  any  nonterminal 
node  of  T corresponding  to  a subdivided  square  sq(u,h,E),  iE|  = 1,2,  will 
be  considered  to  represent  its  center  point  u. 

In  general,  undivided  line  segments  S = sq(u,h,E),  |E|  ■ 1,  do  not  carry 
independent  information.  Thus,  a second  truncation  of  T consists  in  the 
deletion  of  all  terminal  nodes  that  correspond  to  such  line  segments  S. 

With  this,  all  terminal  nodes  of  T represent  undivided,  two-dimensional 
squares  or  points,  and  ell  nonterminal  nodes  correspond  to  points. 


• I 
■} 
:i 


I 
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I 
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Recall  that  in  the  admissible  mesh  defined  by  the  given  tree  T,  a 
[X)int  in  the  open  unit  square  is  irregular  if  it  is  not  a comer 
point  of  all  two-dimensional  squares  incident  with  it  (the  circled  points 
of  Figure  1).  If  conforming  elements  are  used,  as  discussed  in  Section  2.1, 
then  these  irregular  points  carry  no  independent  information.  Obviously, 
an  irregular  point  P can  only  occur  somewhere  on  a divided  line  segment 
S = sq(v,h,E),  |E|  = 1,  which  is  the  side  of  some  undivided  square  sq(u,h,B). 
In  other  words,  any  other  point  on  S must  also  be  irregular.  Thus,  a node 
of  T corresponding  to  an  irregular  point  can  never  have  a descendant  node 
that  represents  a regular  point.  We  may  therefore  truncate  T a third  time 
by  deleting  all  nodes  that  correspond  to  an  irregular  point.  Clearly,  this 
third  truncation  cannot  be  applied  when  nonconforming  elements  are  used. 

By  definition,  all  points  on  the  boundary  aQp  of  the  unit  square 
are  regular.  This  reflects  the  fact  that  on  sQq  other  conditions  have 
to  be  taken  into  consideration.  A side  S of  is  either  the  image  of 
the  intersection  fl  2^  of  two  of  the  subdomains  mentioned  in  Section  2.1 
or  of  a part  of  the  boundary  a2  of  2 where  some  boundary  condition  may 
or  may  not  be  specified.  In  the  first  case  the  points  on  S are  represented 
on  the  subdivision  trees  of  either  one  of  the  subdomains.  Here  it  appears 
to  be  advantageous  not  to  delete  the  corresponding  nodes  from  these  trees 
even  if  they  are  irregular  on  the  union  of  the  subdomains.  The  second  case 
depends  on  the  type  of  boundary  condition.  For  instance,  if  the  solution 
is  specified  on  the  part  of  aQ^  corresponding  to  S then  there  is  cer- 
tainly no  need  to  represent  the  points  of  S on  the  tree.  We  shall  not 
enter  into  the  details  of  the  various  other  possibilities. 


I 


"Aii'iiiLfiaijp 


- 19  - 


At  any  irregular  point  the  solution  is  obtained  by  interpolation. 

Hence  special  consideration  must  be  given  to  the  conputation  of  the  local 
stiffness  matrix  of  an  element  that  has  as  carrier  a square  with  some 
irregular  comers.  This  will  be  discussed  in  Section  3.3.  To  simplify 
this  computation,  a regularity  tag  is  assigned  to  each  node  of  T repre- 
senting a divided  or  undivided  square  Q = sq(u,h,B).  It  consists  of  a 
triple  of  single  bit  numbers  p^  which  indicate  the  regularity 

of  three  of  the  comers  of  Q.  If  the  comer  is  regular,  the  bit  p^  is 
one;  otherwise,  it  is  zero.  The  assignment  of  the  tag  proceeds  recursively 
as  follows: 


(3.1) 


(i)  The  root  of  T has  the  tag  (1,1,1). 

(ii)  If  a nonterminal  node  represents  the  center  u of  a 
divided  square  and  q is  a descendant  node  correspond- 
ing to  the  square  Q,  then  Pq  indicates  the  regularity 
of  the  comer  of  Q opposite  to  u and  p^  that  of 


its  comers  in  the  x^-direction,  i = 1,2,  from  u. 


The  scheme  is  illustrated  in  Figure  5. 

In  Figure  6 we  give  an  example  of 
the  truncated  subdivision  tree  for  the 
mesh  shown  there.  The  nodes  correspond- 
ing to  divided  or  undivided  squares  are 
marked  by  rectangular  boxes.  All  others. 


5; 
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of  course,  carry  no  regularity  tags.  In  all  nodes  the  first  number  is 
the  identifier,  the  pair  of  nimbers  below  it  the  label,  and--where 
applicable- -the  triple  below  that  the  regularity  tag.  It  may  be  noted 
that  the  full  tree  for  the  same  mesh  would  have  more  than  a hundred 
nodes. 


3.2  Subdivision 

Let  q be  a terminal  node  of  T corresponding  to  an  (undivided) 
square  Q = sq(u,h,B).  When  Q is  subdivided  the  following  steps  have  to 
be  taken: 

(1)  New  nodes  have  to  be  created  on  the  tree  corresponding  to  the 
new  squares  created  from  Q and  to  any  points  on  the  sides  of 
Q that  become  regular. 

(2)  Regularity  tags  have  to  be  assigned  to  the  newly  created  nodes 
where  needed;  and  if  some  points  on  the  sides  of  Q have  become 
regular  the  regularity  tags  of  all  nodes  have  to  be  modified  that 
represent  squares  outside  Q incident  with  these  points. 

The  first  part  of  the  resulting  algorithm  creates  the  nodes  for  the 

new  squares  and  assigns  their  Pq  values: 

for  i,j  = tl  do 
begin 

create  nodes  q^j  as  son  of  q; 

XCq^j)  :=  (i.j); 

Pl(qiP  Pjtqij)  :■  0; 
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if  (i.j)  = then  PpC^^j)  ==  1; 

if  (i.j)  = then  PgCqjj)  ==  pQf^)* 

if  (i.j)  = (^i.-^2^  then  Po(qij)  :=  Pi(q); 

if  (i.j)  = (*ii»i2^  then  Pq (q^j)  :=  P2(q); 

end; 

Here  (Xj^.V,)  represents  the  label  of  q. 

The  second  part  of  the  algorithm  needs  the  output  of  the  up-path 
algorithm  of  Section  2.3.  Using  it  we  can  check  whether  the  neighboring 
square  s = s^^  ^ or  s = s^  exists  and  is  divided  in  which  case  the 
center  of  that  side  of  Q becomes  regular  and  a corresponding  node 
t = t^^  0 ^ ~ ^0  +1  created.  Let  v be  the  index  of  the 

coordinate  direction  x , v = 1,2,  and  set  p = 1 if  the  neighboring  square 

is  in  p*"  and  p = 2 if  it  is  in  p . Then  the  down-path  for  a particular 

s or  t starts  from  the  node  where  n = m with  m,  = 2,  m,  = k , 

V = 1,2.  Here  kj,  k2,  of  course,  are  the  indices  (2.20)  obtained  by  the 
uj>-path  algorithm.  We  also  introduce  the  following  label  functions 
of  the  labels  ^2*^2  indices  p,v: 


P V 

1 1 (-x^.il) 

1 2 (±1,-X2) 

2 1 (Xj,±l) 

2 2 (±1,X2) 

Then  the  next  part  of  the  subdivision  algorithm  has  the  generic  form: 


(n,p,a,k)  :=  up-path  (q) 
for  ij.,v  “ 1,2  do 


begin  in  :■  m 


m < n-1  then  s ;=  down-path  (p,n,m,l,2v-3,3-2v) 
if  (ireai-l)  or  (sj< terminal)  then 
begin  for  / = ±1  ^ begin  (i,j)  :=  (Xj ,X2;ti,v) ; 


p (q..)  :=  1 end; 


r down-path(p,n,m,2,v-l,2-v) ; 
create  node  t as  son  of  r; 

X(t)  :=  ((v-1)X^,(2-v)X2); 
if  in  < n-1  then  set  regularity  tags  for  the  squares  in  s 
end; 
end; 


incident  with  t; 


Here  the  condition  m > n-1  indicates  that  the  particular  side  of  Q is 
on  the  boundar>’  dQQ  of  Qq.  Hence  theie  is  no  neighboring  square. 
Nbreover,  since  all  points  on  dQ^  are  treated  here  as  regular  points, 
corres^xinding  nodes  are  always  created  on  the  tree. 

The  statement  requiring  the  setting  of  the  regularity  tags  for  the 
squares  in  s incident  with  t has  the  form  of  the  following  algorithm; 


for  j = tl  do 
begin  for  < = tl  do 

begin  r son  of  s with  label  (-Xj^,-\2;m.,v) ; 

P^(r)  1; 

<diile  (r><terminal)  do 

begin  r :=  son  of  r with  label  (-X^  ; 

;=  ] 


3.  3 tilemcntan-  Stiffness  Matrices 

As  mentioned  in  Section  2.1,  the  appearance  of  irregular  points 
ccmplicates  the  calculation  of  the  elanentary  stiffness  matrices  when 
conforming  elements  are  used.  Iv'e  sketch  here  the  algorithm  for  the  simple 
case  of  conforming,  bilinear,  square  elonents.  It  should  be  evident  how 
tlie  approach  extends  to  more  general  cases. 

Let  Q = sq(u,h,B)  be  a given  square  with  comers  P^,  and 

w.  = w.(x, ,x,;h)  the  usual  bilinear  shape  functions  with  w.(P.)  = 6-, 
i,j  = 1,...,4.  Then  the  desired  solution  y on  Q has  the  form 

4 

(3.2)  y(x^,X2)  = y^w^ (x^,X2;h),  x € Q , 

and--if  all  comers  are  regular--the  elenentary  stiffness  matrix  for  Q is 

(3.3)  ’1^1=  (8(w^,Wj),i,j  = 1,...,4) 

where  B = 8(.,.)  denotes  the  given  bilinear  form.  In  the  case  of  irregu- 
larity some  of  the  values  yj  = y(Pj)  are  given  as  linear  combinations  of 
some  othei  values,  say, 

m 

(3.4)  >"1  “ ^ ^ ^ ^ ‘ 


Then 


n 


s 
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calls  for  the  introduction  of  the  new  functions 

4 


= Wg(x^,X2;h)  = Cj^w'j(x^,X2;h),  I 


- 1 , . . • , 4 , 


and  the  corresponding  elementary  stiffness  matrix  of  Q is  given  by 
(3.5)  M=  (8(w.,w.),  i,j  = l,...,ih)  = C'^MC 

where  C = (c..,  j = 1,...,4,  t = l,...,ft)  is  the  interpolation  matrix  in 
(3.4). 

The  y^  are  the  solution  values  at  certain  regular  points  and  our 
problem  is  to  find  these  points  and  the  interpolation  coefficients  c 
As  an  example,  consider  the  square  4 of  Figure  6.  Here  the  y^  are  values 
at  the  points  5,  13,  21,  22,  and  33  and,  if  the  comers  of  the  square  are 
numbered  counter-clockwise  starting  from  5,  the  interpolation  matrix  has 
the  form 


(3.6) 


(5) 

(13) 

(21) 

(22) 

(33) 

1 

0 

0 

0 

0 \ 1 

0 

1/2 

0 

3/8 

1/8  \ 1 

0 

1 

0 

0 

1/2 

0 

1/4 

1/4 

0 / 

Basically  the  desired  algorithm  is  a modified  version  of  the  down-path 
procedure  of  Section  2.3.  Let  (n,p,a,k)  be  the  output  of  the  up-path  algo- 
rithm starting  from  q.  If  (Xj^,\2)  is  the  label  of  q,  then  the  comer 


I 


P , , (in  the  notation  of  Figure  3)  must  be  regular,  and  the  other 

2 

potentially  irregular  comers  of  Q are  found  on  the  down-paths  from 
p[k^l»  “ 1»2,  with  label  sequences  given  by  (2.21).  These  three  label 
sequences  have  the  general  form 

(3.7)  ie'^  ^ a[<2:iCj-l]}  • a[iCj^:n] 

where,  for  instance,  in  the  case  of  P , . we  have  v = 1,  k,  = k,, 

= 2. 

I, 

Let  R be  one  of  these  comers  and  suppose  that  R is  irregular. 
Tlien  it  must  be  on  a line  S created  by  the  subdivision  of  the  square 
Qf  corresixjnding  to  p[k,].  Then  any  regular  point  on  S is  on  the 
subtree  rooted  at  p[Kj]  with  labels  that  have  a zero  v-th  component. 
Figure  7 shows  a typical  situation. 


pKI  sp* 
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The  branches  deleted  from  the  full  subdivision  tree  are  marked  by- 
dashed  lines.  Only  the  nonzero  comixinent  of  the  label  is  shown  at  each 
node.  Briefly,  we  wish  to  find  the  nodes  and  p ^ on  the  truncated 

A A 

tree  which  correspond  to  the  nearest  regular  points  and  P_j 

bracketing  P. 

The  following  algorithm  for  finding  these  nodes  (if  they  exist)  should 
be  self-explanatory.  It  starts  from  the  node  p[ij],  where  initially 
ij  * and  proceeds  along  the  prescribed  down-path.  If  the  requested 
endpoint  r does  not  exist,  then  the  nodes  returned  with  the 

relative  distances  respectively,  between  the  corresponding  points  and 

Cne  of  these  nodes  may  not  exist- -as  for  the  point  R'  in  Figure  7-- in 
which  case  the  particular  output  node  is  nil. 

Algorithm:  Down- path  interval 
Input  (p,a,ij,i2,v) 

for  j ■ ±1  ^ begin  p^.  :=  nil;  d^  = 0 end 
j - ij;  r :*  p[ij] 
while  Ci>i2)  and  (rj^nil)  do 
begin  i :=  i-1 
n :*  a[i]^ 

A 

p :*  r 

2- 

r :*  son  of  r with  label  a[i] 

end 

^ r - nil  then 
begin  y :■  1/2;  d 0 

while  j > i2  ^ begin  d :■  d +Yali]^;  j i-1;  Y end 

" abs(d);  “ 1 - d,^ 
end 

Output  (r, 


1 
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The  distance  calculation  uses  the  fact  that  by  (2.13) 
dist(R,P  ) -If  ci[j]  2^'^*^ 

where  j is  the  last  index  used  in  the  while  loop.  Since  R must  be 
in  one  half  of  the  interval  centered  at  P we  have 

V 

^ ^ N ^2+i-n 

j-ij 

There  are  t.hree  cases  for  the  output  of  this  algorithm,  namely, 

(i)  r ^ nil,  (ii)  r - nil,  and  no  p^  is  nil,  (iii)  r » nil,  and  one  p^  is 
nil.  In  the  first  case  the  corresponding  row  of  the  interpolation  matrix 
has  a one  in  the  column  corresponding  to  r and  zeroes  elsewhere.  In  the 
second  case,  there  are  the  values  d.j^  in  the  columns  for  p^j  and  zeroes 
elsewhere.  In  the  third  case,  we  have  the  situation  of  the  point  R'  in 

A 

Figure  7;  that  is,  in  the  search  for  one  of  the  points  P^  we  reached  one 
of  the  endpoints  of  S.  Then  we  have  to  apply  our  algorithm  to  the  line 
perpendicular  to  S centered  at  that  endpoint.  The  following  algorithm 
covers  these  three  cases.  The  output  consists  of  the  number  m of  nonzero 
elements  in  the  particular  row  of  C.  These  elements  are  identified  by 

pairs  (r[j],c[j]),  j - 1 vrtiere  c[j]  is  the  element  of  C in  that 

row  and  the  colunn  identified  by  the  node  r[j]. 


aist  - z , 


and  hence 


dist(R,P^) 

dist(P^j,P.j) 


1 

7 


29  > 


Algorithm:  Interpolation  Row 
Input  (p.a.Kj.ic^.v) 
m :«  1;  r[l’  • nil;  c[l]  :»  1 
Vq  • v;  i^:- 
v^ile  r[m]  ■ nil  do 

begin  (r^.p^^.d^j)  down-path  interval  (p,a,ij,i2»''o^ 

^ r nil  then  r[m]  :»  rg 
else  begin 

if  P j^  >=  nil  then  j :*  1 else  j :*  -1; 
r[m]  :=  Pj,  r[m+l]  :=  p_^; 
c[m+l]  :■  c[m]d_j;  c[m]  :*  c[m]dj; 
ra  :*  m+1 

r[m]  * nil  then 

begin  i2  :=  i^;  :=  S-v^;  m :=  a[ij]^  ; ij  = ij+2; 

while  a[i--l]  = m do  i,  :=  i,+l;  ® 

end 
end 
end 

Output  (m,r,c) 

In  the  case  of  r[m]  ® nil  we  must  continue  the  search  on  the  line 
perpendicular  to  the  current  search  direction.  This  is  reflected  by 
the  statement  vq  :=  3 - Vq.  To  find  the  needed  endpoint  of  S is  equi- 
valent to  determining  the  centerpoint  T of  one  of  the  sides  of  Qg  (see 
Figure  3).  Thus,  we  need  to  construct  on  the  up-path  from  correspond- 
ing to  Qg  the  node  from  where  the  new  line  containing  the  desired  endpoint 
is  branching  off.  This  is  the  content  of  the  loop  * ’while  a[ij^'l]^  ■ m do”. 

With  these  algorithms  the  interpolation  matrix  C of  (3.4)  is  constructed 
row  by  row.  The  calculation  of  the  transformed  elementary  stiffness  matrix 

A 

M of  (3.5)  is  then  straightforward. 


3.4  Matrix  Assembly  and  Decoinposition 


The  assembled  global  stiffness  matrix 
has  the  bordered  block-diagonal  structure 
shown  in  Figxnre  8.  The  matrices 
correspond  to  the  finite  element  nodes  in 
the  open  subdomains  i > 1,...,N, 
and  the  border  represents  the  points  on 


In  general,  the  size  of  B is  much 

smaller  than  that  of  the  other  diagonal  blocks;  moreover,  B is  often  not 


very  sparse,  especially  when  N is  small. 

■Hie  form  of  the  matrix  calls  for  the  use  of  block  deconposition.  Since 
the  Aj^ -blocks  are  independent  of  each  other,  it  suffices  to  consider  only 
one  block  at  a time.  The  triangular  factorization 


of  one  of  these  blocks  introduces  the  following  partial  deconposition  of 


the  overall  matrix 


»diere 


(3.9)(ii) 
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Accordingly,  for  each  subdomain  52^^  our  procedure  consists  of  the  four 
steps: 

(i)  Generate  the  matrices  and 

(ii)  Generate  the  contributions  to  the  matrix  B and  send  them  to 
a separate  solution  process; 

(3.10)  (iii)  Con^jute  the  decomposition  (3.8)  of  and  the  corresponding 
modified  matrix 

(iv)  Conpute  the  modification  and  send  it  to  the  solver 

for  B. 

When  these  steps  have  been  conpleted  for  all  subdomains,  the  solver  for  B 
can  conplete  the  deconposition  of  that  matrix. 

The  structure  of  the  solution  routine  for  B is  standard  and  needs  no 
discussion.  Instead  we  consider  only  the  inplementation  of  the  steps  (3.10) 
in  the  setting  of  our  particular  storage  structure.  As  mentioned  before, 
the  tree  structure  allows  for  a natural  use  of  the  nested  dissection  tech- 
nique [ 7 ] in  the  execution  of  step  (iii).  It  also  provides  for  an  effi- 
cient segmentation  of  the  block- row  (A^,C^)  into  groups  of  rows  for  storage 
on  secondary  devices.  The  steps  (i)  and  (ii)  are  combined  with  step  (iii); 
that  is,  the  generation  of  the  coefficients  of  A^,C;  (and  B)  from  the 
corresponding  elementary  stiffness  matrices  is  incorporated  into  the  decom- 
position process.  In  other  words,  an  elementary  stiffness  matrix  is 
generated  (see  Section  3.3)  only  when  it  is  needed  in  the  deconposition 
step  (iii). 


I 
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Any  subdomain  is  the  image  of  the  unit  square  Qq  and  the 
mesh  on  is  defined  by  the  mesh  on  Q^.  Let  T be  the  (tnincated) 
subdivision  tree  corresponding  to  this  mesh.  With  any  node  q of  T 
we  associate  the  first  nonteiminal  node  r^  with  regularity  tag  (1,1,1) 
on  the  path  from  q to  the  root.  For  any  nonterminal  node  q of  T with 
P-tag  (1,1,1)  we  may  then  define  the  node  set  q*  * (p  | p node  of  T with 
tp  ■ q).  By  construction,  this  set  specifies  a subtree  of  T which  shall 
be  denoted  by  T(q*).  Moreover,  in  the  usual  manner  we  derive  frcan  T 
the  c(»itracted  tree  T*  with  the  sets  q*  as  nodes.  Figure  9 shows  T* 
and  one  of  the  subtrees  T(q*)  for  the  tree  of  Figure  6. 


Figure  9 , 
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Note  that  in  T the  three  nodes  q_2»q.i.<lQ  corresponding  to  the 
1 2 

squares  Q_2  * sq(e  ♦e  Q ■ sq(0,2,B)»  and  Qq  have  p-tag  (1,1,1). 

Hence  in  each  case  T*  begins  with  the  succession  of  the  three  nodes 
‘^-2’^-l’^0  •’Wltiple  branching  sets  in.  Of  course, 

r(qf^)  contain  the  nodes  corresponding  to  the  faces  of  Qq. 

Let  N*  be  the  node  set  of  T*  and  N*  ■ N*  ~ 
q*  € N*  let  P(q*)  and  E(q*)  be  the  sets  of  nodes  of  T(q*)  vdiich  in 
the  mesh  correspond  to  points  or  undivided  squares,  respectively.  Because 
of  q ( P(q*)  the  set  P(q)  is  never  empty;  on  the  other  hand,  E(q*)  may 
well  be  empty  as,  for  instance,  for  q*  * q*2.q*i*  The  nodes  of  the  union 
U(P(q*)|q*  € correspond  exactly  to  the  rows  and  columns  of  the  matrix 
block  A.  for  our  particular  subdcmiain  2.  and  those  of  P(q%)  U P(q*i) 
to  those  colunns  of  vrfiich  are  related  to  2^. 

In  order  to  define  the  pivoting  sequence  on  reflecting  nested 
dissection,  we  introduce  botton-up,  left-to-right  orderings  (end-order 
traversals,  see  [8,  p.  334])  on  all  nodes  of  T*  and  on  the  nodes  of  T(q*) 
belonging  to  P(q*): 


(3.11) 


n*:  Ng  - {l,2....,|Ng|} 

qj.*  P(q*)  - {1.2,...,|P(q*)|},  V q*  € N* 


In  the  exanple  of  Figure  10,  q*  and  ^^2  8lve  the  orderings 


18»,20*,22*,23»,27*  and  17,19,21,22  . 
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Note  that  for  the  root  r of  T(q*)  we  always  have  ■HqoCi’)  “ 

The  pivotal  sequence  for  A.  is  now  defined  by  the  lexicographic 
ordering  of  the  pairs 

^ICP)  “ (■n*(p),i1q*(p)),  V p ( P(q*),  q*  ( N*  . 

We  segment  the  rows  of  the  block  row  (A^,Cj^)  into  groups  of  rows 
with  the  same  q*-value 


A(qJ)  C(qJ) 


(A^.C.) 


n 


iNgl 


A(q;;)  CCq;;) 

For  clarity  the  index  i was  siqjpressed  on  the  right  side. 

For  any  q*  C Ng  the  corresponding  node  q^  of  T represents  an 
open  square  with  boundary  lines  q,Tq  and  comers 
as  shown  in  Figure  3 for  Q = Q^.  Let  t^j  q,  t^  and  p^^  be  the 

corresponding  nodes  of  T (if  they  exist).  We  denote  by  q)  and 

^(tfl  ^j)  the  sets  of  nodes  of  the  subtrees  of  T rooted  at  the  indicated 
nodes.  If  the  particular  node  does  not  exist,  the  corresponding  set  is 
empty.  All  (potentially)  nonzero  entries  in  the  upper  triangular  part 
of  the  block  row  (A(q|^)  ,C(qJ^))  then  consist  of  a |P(q^)  | -dimensional 
upper  triangular  matrix  and  off-diagonal  blocks  as  shown  in  Figure  10. 


\ • 2 ^9 

Here  label  of  q^,  and  the  nodes  p ,p  and  p ^ are 

defined  by  the  up-path  algorithm  for  (see  Figure  4) . The  node 
Px^,X2  belongs  to  the  set  PCCqJJ)*)  for  k «=  minCk^.k^.  Under  the 
ordering  (3.11)  the  nodes  of  any  set  ^(t)  / D are  always  consecutively 
nunbered  from  some  index  to  » /q-1  + |tS(t)|.  Here  4q  is  carried 
by  the  left -most,  lowest  descendant  of  t and  by  t itself.  This 
allows  for  a smple  determination  of  the  size  of  the  particular  matrix 
block. 

The  algorithm  for  assembling  and  decomposing  the  block  row  (Aj^,C^) 
now  has  the  following  schematic  form: 

1.  Qear  save  file. 

2.  Establish  the  nodes  of  T*  and  their  ordering  n*. 

3.  Fbr  each  q*  € A/g  establish  the  nunbering  ^^^(p)  of  the  nodes 
of  P(q*). 
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4.  For  q*  € N*  in  order  of  -g*  do 

4.1  Establish  storage  map  for  the  block  row  of  P(q*)  as 
shown  in  Figure  10  with  sizes  deteimined  as  indicated 
above  and  clear  the  matrix  area. 

4.2  q*  is  not  terminal  in  T*  then  read  previously 
generated  updates  for  the  coefficients  of  the  row  P(a*) 
and  delete  them  from  the  save  file. 

4.3  For  all  p € E(q*)  ^ 

4.3.1  generate  the  elementary  stiffness  matrix  M for 
p (see  Section  3.3) 

4.3.2  update  the  row  P(q*)  with  those  coefficients  of 
M for  which  at  least  one  of  the  tivo  index  nodes 
belongs  to  P(q*) 

4.3.3  add  the  coefficients  of  Kl  to  the  save  file  for 

2 

which  ^th  index  nodes  are  in  P((p  )*)  U P((p  )*) 

U P((p  2)*) 

4.4  Decompose  matrix  row  of  P(q*)  and  during  the  decomposi- 
tion add  updates  to  the  save  file  corresponding  to 
coef|icients  in  the  rows  of  P((p^)*),  P((p  ^)*),  and 
P((P  ^)*). 

4.5  Output  decomposed  row  for  P(q*). 


In  the  save  file  (in  secondary  memory)  we  retain  triplets  consisting  of  a 
value  and  two  index  nodes  which  represent  needed  updates  for  matrix  rows  to 
be  decomposed  later.  At  the  end  the  save  file  contains  the  updates  for  the 
matrix  B of  (3.9)(ii).  In  step  3 the  nodes  P(q*2)  riu"'" 

bered.  We  assume  here  that  this  numbering  is  done  outside  of  this  algorithm; 
it  depends  on  the  given  domain  Q and  the  txjundary  conditions.  The  storage 
map  for  the  row  of  P(q*)  establishes  the  (one-to-one)  correspouJonce 
between  the  index  nodes  of  the  rows  and  columns  indicated  in  Figure  10  and 
the  customary  indexing  of  the  matrix  storage  used  here. 


f 
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The  above  algorithm  uses  the  save  file  for  the  full  segmentation 
given  in  (3.13).  In  practice,  several  of  the  block  rows  may  be  taken 
together  using  subtrees  of  T*  for  the  control  of  the  loop  of  step  4. 
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4.  Outlook 

The  general  tree  structure  introduced  here  is  not  restricted  to  the 
particular  meshes  under  discussion.  For  example,  an  extension  to  analogous 
meshes  on  the  closure  Qq  of  the  d-dimensional  open  unit  cube 

Qq  = {X  € I 0 < X.  < 1,  i = 1 d},  d2  1 

is  rather  straightforward.  Here  we  consider  subsets  E c B = {e^,...,e^} 
of  the  set  of  standard  basis  vectors  of  and  define  the  |E| -dimensional 
open  square  Q with  center  u and  side  length  h > 0 by 

Q = sq(u,h,E)  = {x  € r'^  ( ||x-u||„  < Y h,  x^e^  = u^,  V e^  < E}  . 

If  now  the  definition  (2.4)  of  the  "label"  set  A is  changed  to 

A = {a  € R^  I a.  = -1,0, +1,  i = 1 d}, 

then  with  the  same  notations  (2.5), (2.6)  as  before  the  formulas  (2.7)  and 
(2.8)  for  the  faces  of  Q and  for  the  subdivision  operator  a,  respectively, 
remain  exactly  the  same.  With  this,  the  tree  can  be  defined  as  before  and 
the  labelling  rule  carries  over  verbatim.  Thus  also  the  up-path  and  down- 
path  algorithms  extend  to  this  case. 

With  this  generalization  it  becomes  possible  to  combine  trees  for 
meshes  of  different  dimensionality.  For  the  practical  in^lementation  of 
finite  element  con^Jutations,  the  trees  may  again  be  truncated.  However,  in 
higher  dimensions  the  irregularity  problem  becomes  more  complex  since  not 
only  points  but  also  faces  may  be  irregular.  Thus,  the  corresponding 
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truncation  of  the  tree  and  the  interpolation  algorithms  need  further 
analysis.  The  matrix  assembly  and  deconposition  algorithm,  however,  does 
not  change  in  principle. 

The  tree  structure  can  also  be  defined  for  meshes  on  equilateral 
triangles  of  the  form  of  Figure  2 (without  the  irregular  lines).  Basically 
we  need  to  change  the  basis  vectors  in  the  set  B and  take  into  account 
that  now  the  center  point  and  side  length  only  define  the  triangle  up  to 
a rotation  of  size  n.  Clearly,  with  this  an  extension  to  meshes  on  higher 
dimensional  simplices  is  also  possible. 
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